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ABSTRACT 

We provide theoretical procedures and practical recipes to simulate 
non-Gaussian correlated, homogeneous random fields with prescribed 
marginal distributions and cross-correlation structure, either in a N- 
dimensional Cartesian space or on the celestial sphere. We illustrate our 
methods using far-infrared maps obtained with the Infrared Space Obser- 
vatory. However, the methodology presented here can be used in other 
astrophysical applications that require modeling correlated features in sky 
maps, for example, the simulation of multifrequency sky maps where back- 
grounds, sources and noise are correlated and can be modeled by random 
fields. 

Subject headings: methods: data analysis - methods: numerical 

1. Introduction 

Random fields are widely used in astrophysics to model realistic scenarios of 
physical processes that depend on random components. For example, they are 
widely used in cosmology to model different types of galactic and extragalactic 
backgrounds (e.g., Martinez & Saar 2001). The physical characteristics of the 
backgrounds are translated into statistical characteristics of the field. The fields are 
usually required to be ergodic so that information about the model can be extracted 
from a single realization of the field. In addition, they are assumed to be isotropic 
(the correlation between two points depends only on the distance that separates 
them) or homogeneous (the correlation depends on the difference of their position 
vectors) to reflect the geometry of the cosmological model. Given the complexity of 
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the cosmological models, characteristics of the fields are usually studied via Monte 
Carlo simulations where model predictions are compared with observations. 

To simulate a Gaussian random field we only need to specify the mean and 
spectrum (or correlation function) but there are non-Gaussian models whose 
predictions also have to be tested against observations. For example, standard 
inflationary models predict Gaussian temperature fluctuations of the cosmic 
microwave background but there are other cosmological models that predict 
non-Gaussian fluctuations (Avelino et al. 1998; Peebles 1999). Homogeneous 
non- Gaussian random fields are more difficult to define and simulate, and since they 
are not uniquely determined by their first two moments, we often have to accept 
only a partial second order description of the fields. 

Vio, Andreani & Wamsteker (2001) (hereafter VAW) presented numerical 
methods for the simulation of homogeneous scalar random fields R(t) with prescribed 
one-dimensional distribution function (marginal distribution) Fn(r) and correlation 
structure. These methods can be used, for example, to generate cosmic microwave 
background maps with a given spectrum and a marginal distribution that allows 
for asymmetry or kurtosis in the pixel temperatures. In this paper we generalize 
these methods to the case of pair-wise correlated random fields defined on the same 
physical space. This will allow us, for example, to simulate backgrounds and source 
fields that are not independent, as in star-forming regions where the already formed 
stars are linked to the surrounding gaseous and dusty environment from which they 
originate. 

We consider random fields defined in iV-dimensional "parameter spaces" where 
the coordinates t = (ti , t 2 , • • • , ijv) m &y correspond to spatial/angular coordinates 
(spatial random fields), time (time processes), a mix of these two (spatio-temporal 
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random fields), or to even more general situations. We also consider the case of 
random fields defined on the sphere, that is, random fields that depend only on the 
direction in the sky. If the multivariate random field is R(t) = {Ri(t), i? 2 (*) ? - • •> 
-Rm(£)}, the goal is to generate a field with components Ri having (possibly different) 
prescribed marginal distributions F^ri), prescribed correlation functions, and 
prescribed pair-wise cross-correlations. 

In theory, a complete description of a random field requires the definition of 
all finite-dimensional joint distributions, but, unless the field is Gaussian, this is a 
terribly difficult task. In practical applications one only considers a second order 
description of the field that specifies the marginals F Ri (rj), the means [ir z = E[Ri(t) } 
and the cross-covariance functions 



where E[-] stands for expected value (ensemble average). 

As mentioned before, it is often possible, even necessary, to adopt some 
simplifying assumptions like isotropy or homogeneity of the field R(t). In the 
multidimensional context, isotropy means that the cross-covariance function 
depends on the length r = \\t\\ of the vector r = t\ — t 2 but not on its direction: 
{.RiRjiti, £2) = £,RiRj (\\ti — £2 ID- In this case the cross-correlation function (normalized 
covariance function) depends only on r 







where 




(3) 
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and a\. is the variance corresponding to the marginal F^ij-i). By definition, 
P/?iRi (0) = lfor* = l,2,...,M. 

Although the isotropic case is of great interest in astronomical applications, 
here we consider multidimensional homogeneous random fields in general, that is, 
random fields whose covariance function depends on r. In this case, the correlation 
function p R (r) is defined as in equation (2) but with r instead of r. 

The rest of the paper is organized as follows. In Section 2 we describe a general 
procedure for generating non-Gaussian random fields by point-wise transformations 
of Gaussian ones. In Section 3 we provide practical recipes for the methods outlined 
in Section 2 and show how they can be used in either iV-dimensional Cartesian 
spaces or on the celestial sphere. Section 4 presents some examples. We show 
simulated maps of the background and localized source field components of a far-IR 
sky map of the Infrared Space Observatory (ISO). In the ISO map, the locations 
and intensities of the sources are correlated with the background field and only a 
multidimensional simulation may reproduce such physical characteristic. 



2. Transforming Gaussian fields to non-Gaussian 

In Section 3.2 we discuss methods for generating homogeneous Gaussian random 
fields with a prescribed correlation structure. The procedures are straightforward, 
at least in principle, while, in contrast, generating homogeneous non-Gaussian fields 
is not a trivial task. Fortunately, a large family of such fields can be obtained via 
pointwise transformations of homogeneous Gaussian fields. These transformations 
are the basis of the methods used in TZ/HFto simulate non-Gaussian scalar random 
fields. In this paper we generalize them to the multivariate case. 
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We propose the following procedure for simulating multidimensional 
homogeneous non-Gaussian fields R(t) with prescribed cross-correlation structure 
p R {r) and one-dimensional marginals F R .(r): (a) Generate a zero- mean, unit- 
variance, homogeneous multivariate Gaussian field X(t) with a fixed cross-correlation 
structure Px( T )', (b) transform X(t) to a non-Gaussian field R(t) via 

R(t)=g[X(t)], (4) 

where g(-) — {gi(-), g 2 (-), . . . ,9m{')}- The problem is to determine a Gaussian field 
X(t) and functions {gi} that transform Px( T ) into Pr( t )- 

As discussed in VAW, given marginals {Fr.(t)} with no atoms (a concentration 
of a finite probability mass at a point), the integral transform 

Ri(t) = gAX.it) ] = F- 1 ! F x [Xi(t) ] }; i = 1, 2, . . . , M, (5) 

where Fx denotes the standard Gaussian distribution function and F^. 1 is the inverse 
distribution function of Ri(t), provides a random field with marginals {F^^r)}. 
Since </j(-) = F^{F X (-)} is a strictly monotonic function, transformation (5) is 
invertible and pR iRj (r) can be calculated through 

PRiRjir) = / / [gi{xi) - iir.] [gj(x 2 ) - /^] <t>{x u x 2 ; pxiX^r)) dx x dx 2 , 

^Ri&Rj J —oc J — oo 

(6) 

where, X\ = x(t), x 2 = x(t + r), and 

/xx 1 ( + xl - 2p x . x (t) x ± x 2 
PXlXl M ) = 2w[1 _ AixM)r , 2 exp ^ 2|1 _ pU(T)| 

(7) 

For example, Figure 1 shows the relationship between p XiXj and Pr^ for some 
standard marginal distributions. 

Equation (6) can be used to solve for the cross-correlations of a Gaussian field 
that transform to those of the non-Gaussian field. Ogorodnikov & Prigarin (1996) 
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show that pRiRj (t) is a monotonic increasing function of pXiX 3 (t) and that (6) does 
not present invertibility problems. However, although via this inversion it is always 
possible to map a value of pXiXjij) to its corresponding value Pr^^t), there is 
no guarantee that the resulting pXiX 3 (t) is in fact a non-negative definite function 
and hence a valid cross-correlation function. This condition must be checked (see 
below). There are also restrictions on the distributions Fr. that can be obtained. 
For example, as already shown in VAW, the values of pR i R i that can be inverted are 
those in the interval pR i R j £ [Pi^L where p\ and p\ correspond, respectively, to 
the values of —1 and 1 of PXiXj for a fixed function g. This is clearly seen in Figure 
1 where relationships between pR i R i and pXiX 3 for some classic distribution functions 
are shown. 

To conclude this section we note that the fields Ri(t) can be forced to have a 
chosen mean pi and variance of by applying the transformation 

R t (t)=~p Ri + R ^-^ ~a Ri . (8) 

without modifying the cross-correlation structure. 

3. Practical recipes 

To simulate a non-Gaussian field R(t) via equations (5) and (6) we need to (a) 
determine the appropriate cross-correlation functions px^j (t) of the Gaussian fields 
Xi(t) to be transformed, (b) check that Pxi T ) is a non- negative definite function, 
and (c) simulate the Gaussian fields X^t). 
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3.1. Solving the inverse problem 

In principle, pXiXj{r) can be obtained by inverting equation (6) but this 
procedure is computationally expensive. Instead, since the relationship between 
PXiXjir) an d pn^ir) is generally smooth, it is sufficient to choose a set of values in 
the interval [—1, 1] for pXiXjir) and to solve for their corresponding values Pr^r^t) 
via the numerical integration of equation (6). The complete relationship between the 
two correlation functions may be obtained by spline interpolation. The interpolated 
function is then used to determine pXiXjir) as a function of r. 

The function P x {t) has to be non-negative definite to correspond to a 
homogeneous field. This means that every pxiXjir) is a non-negative definite 
function and that the matrix p x {r) is non-negative definite for each r. To check 
this conditions just note that the Fourier transform of pxiXjif), which provides the 
cross-spectrum between the fields X^t) and Xj(t), must be a non- negative function. 
It is therefore sufficient to check that the iV-fold Fourier transform of pXiXjij) 
does not take negative values for any pair In addition, for any fixed r, p x {r) 
has to be a non-negative definite matrix and therefore its eigenvalues have to be 
non-negative. This is automatically checked within the field simulation procedure 
described in the sections below. 

3.2. Simulating Gaussian multivariate random fields in R N 

To simulate a zero- mean homogeneous multivariate Gaussian field X{t) of 
unit-variance components Xi(t) and prescribed cross-correlation functions, we 
use the spectral representation that decomposes the field as a superposition of 
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uncorrelated sinusoids of different frequencies and random amplitudes 

-+oo f+OC 



/-(-oo />H-oo r- 
/ / 
-oo J— oo J — 



dZi(k), 



1,2, 



M; 



(9) 



oo 



N-fold 

where i = \J— 1, k ■ t is the inner product between the wave number and position 
vectors, and dZi(k) are random orthogonal increments satisfying (Priestley 1981) 



E[dZi(k)] = 0, 
E[dZi(k) dZ*(k')\ = 0, k^k', i,j = l, 
E[dZi(k) dZ*(k)} = Sij(k) dk, 



M, 



(10) 



where " * " indicates the complex conjugate operator, dk = dk\ ■ dk,2 • ■ ■ ■ ■ dkw, and 
Sij(k) is the cross-spectrum of the fields X^t) and Xj(t). If the field is Gaussian, 
then dZi(k) is an M-dimensional zero-mean complex Gaussian random variable with 
correlation matrix determined by the cross-spectrum 

/ S 11 (k) S 1M (k) ^ 

s x (k)= ; •. ; , (ii) 

\ S M i(k) S M M{k) ) 
for any fixed k. Increments corresponding to different fc's are independent. 

A similar spectral representations for homogeneous random fields defined on a 
discrete parameter space is (e.g., Ruan & McLaughlin 1998) 



Mt) =££•••£ e^'A^fc); i = 1, 2, . . . , M, 

N v ' 

N-fold 

(the sums are understood to be over all the discrete wave numbers) where 

E[AZ t (k)} = 0, 
E[AZ t (k) AZ*(k')] = 0, k^k', i,j = l,...,M, 
ElAZ^k) AZ*(k)] = Sij(k) Ak, 



(12) 



(13) 
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Ak = Aki ■ Ak 2 ■ . . . ■ Ak N . The discrete representation (12) approaches the 
continuous one (9) as |Afc| approaches zero. As in the one-dimensional case, the 
cross-spectrum Sx(k) is related to p x ( T ) through the (iV-fold) discrete Fourier 
transform. Therefore, under regularity conditions on p x (t), Sx(k) is Hermitian 
and non-negative definite for each value of k, and can be factored as 

M 

S x (k) = H(k)H*(k) = Huim^k). (14) 

i=i 

We can also rewrite AZj(fc) in the form 

M 

AZ^k) = E^(*0 e^^lAfcl 1 / 2 , (15) 
i=i 

where £(k) are independent Rayleigh random deviates and 6i(k) are independent 
phase angles uniformly distributed on the interval [0, 2tt}. We obtain a random phase 
representation of the field by inserting this equation into expression (12) 

M 

Ut) = EE - • -EE e ^ fc ^( fc ) £W e^ (fe) |Afc|^. (16) 

N v '1=1 

N-fold 

Equation (16) can be rewritten in the more revealing form 



Xi(t) = IDFT 



A? 



M 



Y,Hu{k) i{k) 
1=1 



■ |Afc| 1/2 , (17) 



where "IDFTjv[-]" stands for iV-fold inverse discrete Fourier transform. This leads 
to the following procedure for the numerical generation of a discrete field X(t) 
characterized by a prescribed Px( T ) : 

1. Determine the cross-spectra Sx(k) by applying the iV-fold discrete Fourier 
transform to Pxi T )- 

2. Factor Sx(k) as indicated in equation (14). 
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3. Generate a set of M random Fourier increments from equation (15) using a 
set of M independent random phase angles uniformly distributed over [0, 2ir], 
and a set of M independent Rayleigh deviates. The phase angles and the 
Rayleigh deviates at different wave numbers must be independent. Rayleigh 
deviates can be generated by the transformation = [—2 ln-u(fc)] 1 ^ 2 , 
where u(k) is a uniform random number on the interval [0, 1] (Ogorodnikov & 
Prigarin 1996). 

4. Apply the (iV-fold) inverse Fourier transform of the random Fourier increments 
to obtain a set of M complex (or 2M real) random fields. 

Note that the factorization (14) of S x (k) is not unique, there are a variety of 
methods that can be used, some being more stable than others. One possibility is to 
obtain H(k) via the singular value decomposition square root of Sx(k) 

H(k) = S^ 2 (k) = MJ 1/2 M~\ (18) 

where M is a matrix whose columns are eigenvectors of Sx(k), and J is a diagonal 
matrix containing the corresponding eigenvalues. Alternativeley, one can use the 
Cholesky factorization, which may be more numerically stable (Popescu, Deodatis 
& Prevost 1998). With either factorization it is easy to test if the matrix Sx(k) is 
non- negative definite for each k. Indeed, if Sx(k) does not satisfy this condition the 
algorithm providing the Cholesky factorization does not converge, and the matrix J 
contains one or more negative entries. 

Another point to consider is that, although we are interested in the numerical 
simulation of real valued random fields, the suggested procedure does not impose the 
necessary symmetry restrictions on the Rayleigh deviates and phase angles in (16) to 
avoid complex values of the fields Xi(t). However, if the Rayleigh deviates and phase 
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angles are generated with no symmetry restrictions, then the resulting field {Xi(t)} 
is complex with independent real and imaginary parts, each with cross-correlation 
matrix given by p x (r)/2 (e.g., Ruan & McLaughlin 1998). Thus, the real and the 
imaginary parts of \/2Xi(t) provide two independent real valued realizations of the 
desired random fields. 



3.3. Simulating Gaussian multivariate random fields on a sphere 

A zero mean random field R(r) on the unit sphere is homogeneous if the 
covariance of the field between two different directions, r*i and r 2 , depends only on 
the angular distance that separates them. In this case the correlation function can 
be written as 

a RR 

where cos# = r x ■ r 2 , and a 2 RR = E [i?(rx) 2 ] = E [R(r 2 ) 2 ]. The field has the spectral 
representation 

oo I 

R ( r )=Y, Yl a ^ Y ^(r), (20) 
e=o m=-e 

where {Y^ m } is a basis of spherical harmonics. The coefficients { a£ m } are zero 
mean complex uncorrelated random variables with variance that depends only on £ 

E [a e , m al m ] = a 2 £ , (21) 

and < 00 • Similarly, a zero mean multivariate random field R(r) = 

{Ri(r), R M (r)} is homogeneous if each Ri is a homogeneous field and the 
pair-wise cross-correlations depend only on the angular distance. In this case the 
cross-correlation function is as in (2) but with 6 taking the place of r. 

The spectral representation of the cross-correlation function of R is (Yaglom 
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1986) 

p R (6) = J2^ i S e P e (cose), (22) 

where Pi is the Legendre polynomial of degree £ and { Si } are non-negative definite 
matrices. The corresponding spectral representation of the multivariate field is 

R(r) = Y,Ae, m Ye, m (r), (23) 

£,m 

where Ai >m are M x 1 zero mean complex random vectors with 

E [ A f , m A*^ m ] = S e , (24) 
E [ Ai >m A\ m ] = for I 1 or m^m'. (25) 

Each matrix can be factored as in equation (14): = Hn~H\. As mentioned in 
the previous section, different factorizations lead to spectral estimates with different 
variance characteristics. Any such factorization leads to representations like (16) 
and (17) for Gaussian fields in terms of independent uniform random phases 9 s (£,m) 
and independent Rayleigh deviates £ s (£, m) 

R ^ r ) = E [E(%^ m ) erf,(W ] y ^ r ) ( 26 ) 

l,m s 

= SFT[j2(Hi)ksU^™)e i94£ ' m) }; i = l,...,M, (27) 

s 

where SFT stands for spherical Fourier transform. For discrete implementations of 
the spherical Fourier transform see, for example, Dilts (1985) and Driscoll & Healy 
(1994). 

4. Examples 

We start with a simple numerical simulation. Figure 2 shows a simulation of 
three two-dimensional, isotropic, random fields with a standard Gaussian N(0, 1), a 
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Xi, and a uniform U(0, 1) as marginals, respectively. To illustrate different degrees 
of correlation in the simulation, we use the cross-correlation matrix 
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The procedure used is the one described in Section 3 with a matrix H(k) 
obtained by a Cholesky factorization of Sx(k). The marginal distributions in 
the three maps are very different; the Gaussian is symmetric about zero, the x 2 
is nonnegative and right skewed, the uniform is bounded between zero and one. 
However, they all have the same autocorrelation function and there is a clear 
correlation between every pair of maps; the correlation is stronger between the 
Gaussian and the uniform because of the 0.9 coefficient in p R . 



4.1. Simulating the components of ISO far-IR maps 

VA W presented simulations of ISO far-IR maps but their method was unable to 
integrate the correlation between the field of localized point sources and the diffuse 
background. This integration is needed since in the observed data the location 
and the intensity of the brightest sources are correlated with the brightest regions 
in the background. As a result, the methodology presented in WAV is unable to 
satisfactorily reproduce the localized source features of the original ISO map on 
large scale. 

We now show that the methods presented in the previous sections can account 
for this type of correlation information and efficiently exploit it in generating the 
simulated map. Figure 3 shows simulations of the random background and source 
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components of a far-IR map obtained by the following procedure. Figure 3a is the 
observed map, recorded by the ISOPHOT camera (Lemke, Klaas & Abolins 1996) 
on board the ISO satellite (Kessler et al. 1996) at 175 /zm, with a projected size 
of roughly 40' x 40' (Dole, et al. 2001). This map is used to generate a Gaussian 
background field (Figure 3b) via Fourier phase randomization. If physical models of 
the diffuse background were available, they could also be used to simulate this field. 

It is well known that extreme values of distributions can be used to generate 
point processes. For example, consider Xi,...,X n uncorrelated Gaussian random 
variables. Define N as the number of Xi above a threshold r and p = P[Xi > r]. 
Then N is a Binomial B(n,p) that can be approximated by a Poisson with A = np for 
n large and p small. In our case we need more than a point process for the location 
of the sources, we also need a distribution for their intensity. For this reason we have 
adopted a different approach based on Gaussian mixtures. These distributions have 
been successfully used to model outliers and other types of nonstationary behavior in 
time series (McLachlan & Peel 2000; Le, Martin & Raftery 1996). Recall that the 
distribution of a A;-component Gaussian mixture is of the form F(x) = X^=i a iGi(x), 
where the on > (£V on = 1) define the probabilities of sampling from each of 
the Gaussian distributions Gi (note that a Gaussian mixture with more than one 
component is not Gaussian). 

For the ISO example we have used a mixture of two Gaussian components for 
the marginal distribution of the pixels. The idea is the following. The location of the 
extreme values of this mixture define the point process, the intensities are defined 
by the Gaussian distribution of the first mixture component. The second component 
generates an essentially zero background. For example, Figure 4 shows a realization 
of the Gaussian mixture used for Figure 3c (the intensities in the original ISO map 
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were normalized to unit standard deviation). We also define a cross-correlation 
function to correlate the location and intensities of the sources with the Gaussian 
background in Figure 3b. 

Since we do not have a separate image of the source field, in order to fix the 
parameters of the mixture, we identified point sources by choosing peaks in the 
original ISO image. The mean and standard deviation of the selected sources were 
used to define the Gaussian component with a fraction a equal to the proportion 
of pixels identified as sources. We used a Gaussian centered at zero with a small 
standard deviation for the second mixture component (see Figure 4). This means 
that a proportion 1 — a of the simulated observations are essentially zero. The 
rest of the observations have an intensity spread about the mean of the observed 
sources. In practice the Gaussians are truncated to avoid numerical instabilities in 
the integral transform (5). 

Finally we have to choose a cross-correlation matrix. To illustrate the 
correlation of the sources with the bright regions in the background, we used the 
cross-correlation function 



where p is the auto-correlation function of the original ISO map, (3 G [—1, 1] is 
a constant, and B{r) is the correlation function introduced by the instrument's 
smoothing. The choice of this form for p(r) is based on a lack of sufficient information 
about the source field. In particular, the cross-correlation of the background and 
source fields is assumed to be a multiple of po only for illustrative purposes; there are 
of course applications where this is not true, see for example Mo & White (1996). 
For the same reason, we decided to assume that the correlations in the source field 
were introduced only by the Gaussian beam smoothing of the instrument. Figures 




(29) 
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3c, d show fields of sparse localized sources for f3 = 0.20 and f3 = 0.90, respectively. 
In both cases, the check described in Section 3.1 has confirmed that p(r) is a valid 
cross-correlation matrix for each r. It is clear that the correlation of the locations 
and intensities of the sources with the background map of Figure 3b increases with 
(3. 

The application of this methodology for the realistic simulation of far-IR maps 
recorded by the array cameras on board the Herschel satellite (Pilbratt 2001) is 
under investigation. The combination of both physical models and observations is 
compelling in determining the types of distributions and cross-correlation functions 
that are appropriate for the particular application under investigation. 

5. Summary 

We have described spectral methods for the simulation of homogeneous non- 
Gaussian multivariate random fields with prescribed cross-correlation and marginal 
distributions. The basic idea is to apply pointwise transformations to homogeneous 
Gaussian fields to guarantee the homogeneity of the field. The transformations are 
chosen to obtain the desired marginal distributions and cross-correlation function. 
The proposed methodology works for a wide range of cross-correlation functions and 
marginal distributions of general interest. Among these, some Gaussian mixture 
models may be useful to model fields of localized sources. 

We have illustrated the methods by simulating sky maps with background and 
source components that are physically linked. The same methods can also be applied 
to simulate and/or analyse more complex maps, such as those from all-sky surveys 
of the Planck mission (Mandolesi et al., 1998; Puget et al., 1998). The objective of 
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the Planck Satellite is to observe the Cosmic Microwave Background anisotropics 
which account, however, for only a very small part of the much more conspicuous 
backgrounds, either Galactic and Extragalactic. Sky maps will be produced at 
nine different frequencies but strong correlations are expected among different 
maps because of the wide frequency distribution of the expected emissions. Many 
algorithms have been already tested to recover as much information as possible 
from the observed maps (see e.g. Maino et al., 2001 and references therein). But to 
extract important cosmological information the methods should provide information 
on all the contributing components and should be able to model the correlations 
among them. The methods we have presented here are based on theoretically simple 
assumptions but do require data, or other available physical information, to select 
the correct model parameters (cross-correlation function, marginal distributions, 
Gaussian mixture parameters, etc.). Appropriate modeling of more complex sky 
maps is the goal of future research. 
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Fig. 2. — Realizations of three correlated random fields with prescribed marginal 
distributions and cross-correlation structure given by equation (28). The correlation 
function p°(k) is also shown for reference. 
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Fig. 3. — Panel (a) shows the original ISO far-IR map. The simulated Gaussian 
background map obtained via Fourier phase randomization is shown in (b). Panels 
(c) and (d) show simulated source fields for a = 0.20 and a = 0.90, respectively. 
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Fig. 4. — One-dimensional realization of a sequence of uncorrelated deviates from 
the two-component Gaussian mixture used for Figure 3 (see text). The intensity was 
based on the ISO map normalized to unit standard deviation. 



